Phase Diagrams for Spin-1 Bosons in an Optical Lattice 
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In this paper, the phase diagrams of a polar spin-1 Bose gas in a three-dimensional optical 
lattice with linear and quadratic Zeeman effects both at zero and finite temperatures are obtained 
within mean-field theory. The phase diagrams can be regrouped to two different parameter regimes 
depending on the magnitude of the quadratic Zeeman effect Q. For large Q, only a first-order phase 
transition from the nematic (NM) phase to the fully magnetic (FM) phase is found, while in the case 
of small Q, a first-order phase transition from the nematic phase to the partially magnetic (PM) 
phase , plus a second-order phase transition from the PM phase to the FM phase is obtained. If a 
net magnetization in the system exists, the first-order phase transition causes a coexistence of two 
phases and phase separation: for large Q, NM and FM phases and for small Q, NM and PM phases. 
The phase diagrams in terms of net magnetization are also obtained. 

PACS numbers: 37.10. Jk,03.75.-b,75.25.+z 
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I. INTRODUCTION 

The study of cold atoms in optical lattices has cap- 
tured a lot of recent attention. A primary motivation is 
to study the strongly repulsive (two spin species) Fermi 
Hubbard model in the regime of close to one atom per lat- 
tice site in two dimension, a system which is believed by 
many to capture the most essential physics of the high 
temperature oxide superconductors [1]. Much progress 
has already been made towards this goal, in particular 
the Mott insulating phase in three dimension has al- 
ready been obtained @,0|- However, the expected anti- 
ferromagnetic Neel ordering has not yet been reported, 
perhaps due to the difficulty in cooling fermions. 

On the other hand, there are also substantial in- 
terests in studying Bosons with spins in the Mott in- 
sulating regime in an optical lattice. There have al- 
ready been quite a number of experimental studies on 
spinor Bose-Einstein condensates (without optical lat- 
tice) 0, i, H, 0, H, H- Mott insulating state of Bosons 
with frozen spin degree of freedom has also been achieved 
experimentally. [1Q] Hence, one can be hopeful that we 
can study experimentally Bosons with spin in an opti- 
cal lattice in the Mott regime, where though there is 
no net mass transport possible, the spin degree of free- 
dom is still active. Due to the finite tunneling amplitude 
and hence exchange interaction between bosons on neigh- 
boring sites, one again expect the possibility of study- 
ing quantum magnetism and ordering in these systems. 
Moreover, it can easily be seen that the spin Hamil- 
tonian realized in these systems would be very differ- 
ent from their counterpart in solid state magnetic sys- 
tems. For example, for spin-1 atoms, the Hamiltonian 
coupling neighboring spins S.y is of the form [ll|, [l2j 
J(Sj • Sj) + K(Si ■ Sj) 2 with K of the same order as J. 
This is very different from the usual Heisenberg Hamil- 
tonian J (Si ■ Sj) which well describes electronic spin in- 
teraction in solids. Indeed, a large number of theoretical 
papers have already been devoted to the subject of the 
spin physics in these systems, (see [IE, EE EE EE EE Ell 



and references therein). 

In this paper, we consider spin-1 Bosons in an isotropic 
three-dimensional optical lattice in the Mott regime of 
one particle per site. We are in particular interested in 
the case of anti-ferromagnetic interaction between the 
atoms, as in the case 23 Na. The Hamiltonian [H], E2 
correspond to J < 0, K < with |J| < This spin 
Hamiltonian has already been considered in the literature 
even before the field of cold atoms [lEIIlJ- A general con- 
sensus was that, at low temperatures, the system would 
order in a nematic state which breaks rotational symme- 
try but has no net spin on any site; (The dimer state, 

the ground state in one-dimension [II], EE], is unstable 
towards the nematic state with sufficiently strong cou- 
pling between neighboring chains [IE]). However, there 
are some issues in cold-atom systems which were not con- 
sidered in these works, and we would like to remedy a few 
of these in this paper. One is the existence of finite mag- 
netic fields in realistic experiments. This magnetic field 
produces a "quadratic Zeeman" effect [2l|. which lifts the 
energy degeneracy between two atoms in the rrif = hy- 
pcrffiie sublevel versus one each in ro/ = ±1. The other 
consideration is that, in the time scale of the experiment, 
the net " magnetization" , namely the sum of m/ over all 
the particles, is conserved. This "constant magnetiza- 
tion" constraint was usually ignored in previous stud- 
ies. Since in particular the nematic state itself carries no 
magnetization, it is natural to ask what is the thermody- 
namical state of the system if one is constrained to have 
a finite net magnetization. Besides intrinsic interest, this 
issue may be relevant since a realistic experiment may not 
always have exactly equal numbers of my = ±1 atoms in 
its initial preparation. Lastly, one need to consider finite 
temperatures. The nematic state can now tolerate some 
net magnetization via thermally excited particles, and it 
is of interest to know what this amount would be. 

In a previous paper [221 ]. we have already considered 
the finite temperature thermodynamical properties of the 
nematic state, but without the effect of finite magnetiza- 
tion and quadratic Zeeman field. There we in particular 
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have evaluated the entropy of the system, and showed 
that the nematic state can tolerate a large entropy with- 
out being disordered. Since it is now routine that Bose- 
Einstein condensates be cooled to very low temperatures, 
it should therefore be relatively easy to reach this ne- 
matic state by ramping up an optical lattice from a Bose- 
Einstein condensate. We are therefore particularly hope- 
ful that physics of the mentioned spin Hamiltonian can 
be studied in the cold-atom systems. 

For the reader's convenience, the different phases con- 
cluded in this paper are pictorially shown in Fig. [T] for 
zero temperature and in Fig. [2] for finite temperatures. 
At zero temperature, the phases depend on the ground 
states. For larger magnitude Q of the quadratic Zeeman 
effect (we shall provide the condition how large Q should 
be in the main text) , only two kinds of states appear: the 
nematic (NM) state with zero magnetization per site m 
and the fully magnetic (FM) state with m = 1, as shown 
in Fig. Ufa). In between the two states coexist and are 
spatially separated. For smaller magnitude Q the two 
states remain for m = (NM) and m = 1 (FM), how- 
ever, a new state appears above the magnetization m m i n : 
the partially magnetic (PM) state, as shown in Fig. [T](b). 
This new state breaks the rotational symmetry along z 
axis and has magnetization smaller than 1. If the mag- 
netization is between zero and m m i n , there coexist the 
NM state and the PM state. For finite temperature, the 
phase pictures are slightly changed as shown in Fig. [2] 
The system is not a pure state anymore, but a statis- 
tical mixture of different states. For larger Q, we have 
the NM phase if the net magnetization in the system is 
between and a small value mi, while the FM phase is 
obtained if mi < m < 1. In between, phase separation 
of NM and FM phases is expected. This is shown in Fig. 
Q] (a). On the other hand, if Q is small, the PM phase 
will appear as at zero temperature. The NM phase ap- 
pears with very small magnetization m < m,\. The PM 
phase appears spatially separated from NM above m\ 
and occupies an increasing volume fraction with increas- 
ing magnetization. When m m i n reached, the PM phase 
occupies all the region. Above m,2, the system is in the 
FM phase. 

Our paper is organized as follows. In section |TT] the 
model for a strongly repulsive atom-atom interaction in 
an optical lattice with linear and quadratic Zeeman ef- 
fects is introduced. In section HTTl we provide a mean-field 
treatment to solve the problem. In section HVl the phase 
diagrams are obtained either as a function of magnitude 
of linear Zeeman effect or as a function of the magne- 
tization, both at zero temperature PV Ap and at finite 
temperatures (|IVB[) . In section [V] some additional dis- 
cussions and the conclusion are made. 



II. MODEL 

In this paper, we consider spin-1 Bosons loaded in an 
strong optical lattice under the influence of linear and 




FIG. 1: (color online) Phase pictures for zero temperature. 
| shows the nematic (NM) state, f the fully magnetic (FM) 
state and f the partially magnetic (PM) state, (a) For larger 
Q. m = 07 the NM phase; < m < 1: NM and FM phases 
coexist and are spatially separated, m = 1: the FM state, 
(b) For small Q. m = 0: NM; < m < m min : NM and PM, 
phase separation; m m in < m < 1: PM and m = 1 FM. 



quadratic Zeeman effects. In the case of one atom per po- 
tential well, such systems can be described by the Hamil- 
tonian 



H? + H 



(1) 



<i,3> 



where the two-body Hamiltonian H. L j is related to Bose 
Hubbard model and < i,j > denotes the next-neighbor 
sites. Defining the hopping constant t and the interaction 
strength Us depending on the total spin S = 0,2, the 
on-site repulsion coefficients in Bose Hubbard model, the 
energy of the two-body system can be classified according 
to the total spin and therefore can be written as 

01113,122 



Ha = e P, 



(o) 



(2) 



(2) 



where eo — — e2 = — anc ^ P r °j ec ti° n opera- 

(S) 

tors Py project the pair i,j into a total spin hyperfine 
spin S state. The if,- term results from magnetization 
conservation and the linear Zeeman splitting [2 1\ 



(3) 



and iJ 4 is quadratic Zeeman Hamiltonian 

H? =4Q (n u + ni,i) 



(4) 

with n-f , no and n± representing the number operators 
with S z = 1,0, —I, respectively. The two-body Hamil- 
tonian Hij can also be written in a spin representation 

Hj = ■ Sj) + K(Si ■ S,) 2 + J-K, (5) 
where J = ea/2, K = (2eo + e 2 )/6 . 
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FIG. 2: (color online) Phase pictures for finite temperatures. 
Different symbols are shown. [: NM. f: FM with positive 
magnetization, j: FM with negative magnetization, f : PM 
with positive magnetization. J.: PM with negative magneti- 
zation, (a) For large Q and (b) for small Q. They are quite 
similar to the zero-temperature case as shown in Fig. [T] How- 
ever, there still exist some differences. At finite temperature, 
the system is composed of statistically mixed states, for ex- 
ample, the NM phase consists of mostly NM states and but 
mixed with small amounts of FM states. The PM phase con- 
sists of large amounts of PM states and small amounts of NM 
states, etc.. The second difference is that the NM and FM 
phase have a range of magnetization due to this statistical 
mixture. 



III. MEAN-FIELD TREATMENT 

As we have mentioned in our recent paper 22], in order 
to describe the broken 0(3) symmetry for nematic state, 
one can define a new set of basis, 

I*) = ^(-It> + li», 

\y) = ^=(IT> + li», 

\z) = |0) 



(6) 
(7) 



In this basis, the two-body Hamiltonian can be ex- 
pressed as a sum of zero order and second order polyno- 
mials 



Hij = e Pg ] + e 2 pg ] 
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E 

a.0—x : y,z 



a); a 



e2$>></|, 
{i} 



(8) 



where \I a0 ) = 4^(\a)i\P)j + \f3)i\a)j) for a ^ (3, 
l J o> = \/f(l z )*l 2 )j ~ \\ x )i\ x )j ~ \\v)i\y)j) and \h) = 



jK{\x)i\x)j — \y)i\y)j), and the linear and quadratic Zee- 



man Hamiltonian have the form 

Hi = -i^(\y)u(x\ - \x)u(y\) 



H 



Q 



±Q{\x)u(x\ + \y) u (y\). 



(9) 



(10) 



Without H L and terms, we have seen that 
the density matrix should have the diagonalized form 
^2a=x y z P aa \&)((x\- This obviously remains valid when 
HQ is included. However, we see that H L contains off- 
diagonal terms in the |a:),|j/) representation. Therefore 
the general density matrix of a single site should have 
the form 



P = 



E p° 

.=x,y,z 



\a)(c 



P xv \x)(y\ 



P yx \y){* 



(ii) 



where p aa are real and {p xy )* = p yx due to the hermitic- 
ity of p. In this way p xy and p yx can be chosen purely 
imaginary because together with p xx and p yy the real 
part of p xy and p yx forms a real symmetric matrix and 
therefore can be diagonalized. In other words, one can ro- 
tate the system along z-axis to make p xy — —p yx = —ip" 
with real number p" . 

The principle of mean field theory is to reduce a many- 
body problem to a one-body problem by replacing all in- 
teractions to any one body with an average of effective 
interaction. A mean- field treatment for a spin-1 Bosons 
in a lattice has been done by different authors [lj| • 
For Hamiltonian |T]), the only term which has to be av- 
eraged is the two-body Hamiltonian Hy. The effective 
Hamiltonian to replace be a single site operator 

H°eff = ^[pjHij] = zTrj \ Pj (e P^ + e 2 P t 



>(2)x 

ij j 



(12) 



with the coordinate number z. For a cubic three- 

1° 



dimensional lattice, z = 6. Using Eqs.© and (fTTj) , H° 



can be obtained as 



H° eff =z (Kp aa + ^)\arn 



a.—x,y,z 

-z(2J-K){p* y \x)(y\+p yx \y)(x\} 



(13) 



The total effective Hamiltonian H e f f has to include the 
linear and quadratic Zeeman effect as well 



Heff =H e ff + H 



=H° eff -i\{\y){x\-\x){y\} 



(14) 



+ AQ{\x)(x\ + \y)(y\}. 
Defining a new set of parameters h a @ as 
H*f f = - Yl h aa \a)(a\~h x y\x)(y\~hy x \y)(x\, (15) 



V2 



a—x,y,z 
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and comparing Eq. ([15]) with Eqs. (|T3]) and (| 14[) we obtain 



fl 


= -z[Kp + —) 


— 




h yy 


= -z(Kp yy + ^) 


-4Q 




h zz 


= -z(Kp?> + ^) 




(16) 


h xy 


= z{K - 2J)p xy - 


i\ = -iftll 




h yx 


= -h xy . 








ft 11 = z[K-2J)p 


+ A 


(17) 



where 



and all other components are zero. h a @ can be one-to- 
one mapped into p a @ and therefore we can use h a @ as 
parameters to find self-consistent equations for the mean- 
field theory. 

To find the self-consistent equations we first rewrite 
H e ff in a matrix representation in (\x), \y), \z)) T basis. 
H e ff therefore has the form 



H eff = 



^ hid,, 



(18) 



where 





" 1 
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1 
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1 
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^3 = 





-1 





T z = 






























1 



and 



foxx _j_ foyy 

ho = , hi = 0, h 2 



h xx — ftyy 
hz =- , h z 



(19) 



(20) 



h* 



For the convenience of latter use, we can also define pi 
and p z in the same way 



Po 

P3 



p xx + p 



,</</ 



-, Pi = 0, p 2 = p\ 



p xx - p yy 



(21) 



Pz = P 



The one-body density matrix of canonical ensemble is 
defined as 



exp(-/3ff e// ) 
Trexp (-/3H eff ) 



(22) 



where (3 = l/k B T. Inserting Eqs.CEH]), (0 and $ZUj) into 
Eq. (|22|) and after some algebra (see Appendix A), p 
reads 



P0°0 + P2&2 + PZ&Z + PzT z 



(23) 



where 
Po = 

P2 = 



cosh /3 ft 



P2 



I12 sinh (3h 



e &h zo + 2cosh/3/i' 

hz sinh j3h 
h{eP h -° + 2cosh/?ft)' 



ft(e 



2 cosh /3ft) 



2 cosh /3ft) 

(24) 



with the definitions: ft, = ^h^ + h^ and h zo = h z — ho- 

Comparing Eq. (Tl7o|) with Eq. (j2"3"|) , one can obtain 
three self-consistent equations through the definition of 
h ZOl h2 and hz in Eq. (|20[) . The first equation can be 
obtained by the relation h zo = — zK(p z — po)+AQ, which 
leads to 



z\K\ 



cosh /3 ft 



O f3h zo 



2 cosh /3 ft. 



■4Q. 



(25) 



/i 2 (= z(i^T — 2J)p 2 + A) accounting for the off-diagonal 
term in the effective Hamiltonian gives the second equa- 
tion 



h 2 = z{K-2J) 



ft.2 sinh f3h 



h( e P h *° +2 cosh /3ft) 



+ A. 



(26) 



The third equation can be found by the relation: ft 3 
z\K\pz, which gives the form 



hzz\K\ 



sinh [3h 



ft(e^° + 2cosh/3ft)' 



(27) 



Therefore there are two situations: if hz is nonzero, then 
Eq. (|2T|) can be reduced to 



sinh /3 ft 



z\K\ (eP h »° +2 cosh /3ft)' 
Inserting Eq. (|2"8)) into Eq. (|2"6"|) , ft 2 is a constant 

\K\X 



2A 



(28) 



(29) 



with the definition: A = J — K > 0. In the case that 
hz = 0, ft = ft-2 and Eq. (|2"6")) is also reduced to a two 
parameter equation 



ft 2 = z(K - 2 J) 



ft-2 sinh (3h,2 



ft(e 



ph zo 



2cosh/3ft 2 ) 



+ A. 



(30) 



In either case we have reduced the mean-field problem to 
two self-consistent equations. 

These self-consistent equations may have many solu- 
tions, however, only the one which has the lowest free en- 
ergy describes the equilibrium state of the system. There- 
fore we should find the free energy. The free energy can 
be calculated by the relation 



F — Ej, 



E e 



TS 



(31) 



where the internal energy is given by the two-body in- 
teractions Ei n t = \ Tr pH^ff and the external energy 
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FIG. 3: (color online)Energy difference of the NM, FM and 
PM states with respect to the NM state vs. A for 4Q > zA/2. 
The black line shows the reference state: NM, the blue line 
shows the FM state and red-dashed line shows the PM state. 
A n j is the point of a first order phase transition. In the inset 
we show the phase diagram in terms of A. 



is given by the Zeeman fields E ext = Tr p{H L + H®} . 
After some algebra (see Appendix B) the free energy is 
obtained as follows: 



zK 



■2zA 



4e /3,l "cosh/3/i + 2 " 
~ (e^-o + 2cosh/3/i) 2 _ 

h\ sinh 2 (3h ze 2 
hpf^hTo + 2 cosh /%) 2 + ~'' 



(32) 



E„ 



4Q 
2A 



a /3h, c 



e ph zo + 2 C o S h (3h 
h 2 sinh j3h 



(33) 



h(e 



2 cosh (3h) 



4Q 



and 



-TS 



h zn e? h * 



2h sinh (3h 



eP h *° + 2 cosh 0h e$ h * 
1 



2 cosh (3h 



(34) 



\TL{eP h *° + 2cosh/3/i). 



IV. PHASE DIAGRAMS 

In the following, we discuss different phases at zero 
temperature and at finite temperatures. 



A. Zero Temperature 

At zero temperature, p e ff is dominated by the smallest 
eigenvalue of H e f f , which can be easily found by diago- 
nalizing Eq. (| 14|) . H e ff can be rewritten as follows 



H e ff = -h l - h 2 a 2 - h 3 a 3 - h zo r z 



(35) 



with identity matrix 1. Obviously hot is a constant ma- 
trix, therefore it can be ignored. We define a new Hamil- 
tonian 



TTll 

H eff 



-h%cr%- h 3 <r 3 - h zo r z 



(36) 



One has three eigenvalues for Et^ff'- eigenvalue — h zo 
correspond to eigenvector [0, 0, l] T (|z)) and eigenvalues 
=fh correspond to eigenvectors: [u±, v±, 0] T , i.e. (u±\x) + 
v ±\y))- u-\x) + v-\y) is irrelevant at T = because its 
eigenvalue h is positive. Therefore if h zo > h, the system 
is in the pure nematic state \z), otherwise the system is 
in the (u + |x) + v + \y}) state. This state can be either a 
patially magnetic (PM) state or a fully magnetic (FM) 
state depending on the parameters A and q. We will 
discuss the details later. 

In the case h zo > h, the free energy can be calcu- 
lated by using Eqs. (|3"Tj) - (|34| . We can see that at zero 
temperature the equations show the competition between 
h zo and h. After some algebra, F can be rewritten as a 
function of coshj3h/e h *° and of sinh f3h/e hz ° . These two 
terms disappear at zero temperature. Therefore 



F 



zK 



ze 2 
T 



(37) 



zK 
2 



-^f 2 - = z(e +2e2)/6 is independent 



We see that E z 
of A and Q. 

On the contrary, if ft, > h zo , the eigenvector of El^j 
can be solved by the equation below 

(h 2 <J2 + h 3 a 3 )(u+\x) + v+\y)) = h(u+\x) + v+\y}). (38) 



This yields 



ih 2 



1 



1 h 3 i 

~72 [ + T ] ^ V+ - V2h(l+<f)h 



(39) 



In a similar way, by using Eq. (|31|) - (|34[) , the free energy 
defined as E + has the form 



F 



zK zeo 



zA fhi 
2 \h 



(40) 



We can define h 2 /h = cos 9 and h 3 /h = sin 9 since h 2 
h\ + h 3 . E + then takes the form 



E 4 



zK 
~2~ 



ze 2 



zA 



cos 2 9 - \cos9 + 4Q. (41) 



There are two minima for E + : either 
sin 9q = 



or 



~zA' 



(42) 



(43) 



The second solution has a constraint: A < zA, otherwise 
there is no solution due to the fact that cos# can not 



be larger than 1. These two saddle points can be also 
obtained by the self-consistent equations. In the case 
/13 = 0, this indicates directly that sin6> = 0. This yields 
8q = and then u = and v = according to Eq. 

([55)) . Therefore the ground state reads 



\*) = -j={\x)+i\y)) = -\]). 



(44) 



Therefore we obtain a fully magnetic (FM) state . On 
the contrary, if hi 7^ at T = 0, Eq. ([28| is reduced to 
the form 



z\K\ 



(45) 



Together with Eq. ([25]). we obtain Eq. 03]). The eigen- 
state of this solution is 



I*) = u\x) + v\ 



where u, v(^= ±1) are given by Eq 
state into spin basis, we obtain a state 



I*) 



a| T>+7l I) 



(46) 

Transforming the 
(47) 



where 



and 



u + iv 1 ho , 



(48) 



7 : 



V2 



1 
71 



1 



(49) 



by using Eq. ([59"]) . We call this a partially magnetic (PM) 
state. We note that /13 7^ implies that x and y axes are 
no longer equivalent and the rotational symmetry about 
the z-axis is spontaneously broken in this PM state. 

We can summarize that we have three phases: nematic 
state (NM) \z), FM state | T) and PM state a\ T) +7| ! 
). To see which state is preferred we have to calculate 
the free energy for these three states. Define the energy 
difference first : AE = E + — E z . This yields 



AE(9) = ^(cosfl) 2 



Acos6> + 4Q. 



(50) 



Therefore a FM state has the energy difference to a NM 
state 



AE(d ) 



zA 
~~2 



4Q- A, 



while a PM state has the energy difference 

A 2 



AE{6 1 )=AQ- 



2zA 



(51) 



(52) 



with the constraint: A < zA. A state is favored over the 
NM states only if AE < 0. Therefore a PM state can be 
a ground state if there exists a critical lambda X n p 



A,, 



\/8QzA 



(53) 



AE 




FIG. 4: (color online) Energy difference of the NM, FM and 
PM states with respect to the NM state vs. A for < 4Q < 
zA/2. The black line shows the reference state: NM, the blue 
line shows the FM state and red-dashed line shows the PM 
state. X n ,p is the point that a first order phase transition 
occurs from NM to PM and X p j is the second-order phase- 
transition point from PM to FM . In the inset we show the 
phase diagram in terms of A. 



where X n>p < zA.That means a PM state can be a ground 
state only with the condition 



zA 

4Q<-. 



(54) 



This separates the whole parameter space into two 
regimes: a regime with PM states and a regime with- 
out. 

(a) 4Q > ^y'- In this regime, there exist only two 
states: NM and FM. Fig. [3J shows AE of different states. 
Since we subtract the energy of the nematic state in the 
definition of AE, we can define AE = for the nematic 
state, as the black line shown in Fig. [3J The blue line 
decreasing linearly shows AE(0q), the energy difference 
for FM ([5T]) . AE(8q) becomes negative if A < X n j, where 



A 



n,f 



4Q 



zA 
~2~' 



(55) 



The system undergoes a first-order phase transition from 
nematic states to fully magnetic states while A passing 
A n /. This picture is also drawn in Fig. [3] The rea- 
son why the phase transition is first-order is that the 
magnetization jumps from zero for nematic states to 
one for ferromagnetic states. We note that from Eq.Q 
dAE/dX = — (n.| — nj,), hence the slope of AE versus 
A is proportional to the magnetization. In order to see 
that the PM state does not appear in this regime, we also 
draw AE(6 1 ) as the red line in Fig. [3] AE(9x) is always 
positive till the end point A = zA. Therefore PM never 
appears in this regime. 

In experiments the magnetization is constant in time, 
therefore it is important to have a phase diagram with 
magnetization as a parameter. Supposed that average 
magnetization per site is m, the system is purely NM only 
if m = 0, while it is purely FM only if m = 1. In between 
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(a) 



NM 



NM and FM , PS 



FM 



NM and PM , P.S. 



PM 



1 

FM 



FIG. 5: (color online) Phase diagrams for zero Temperature 
in terms of average net magnetization m. (a) 4Q > zA/2 
(b)0 < 4Q < zA/2. P.S. means phase separation. This figure 
corresponds to Fig. [T] 



the nematic state and the PM state spatially separated. 
Supposed that the fraction in nematic state is defined as 

x, we obtain x = 1 — m y/§0- In the regime: m m in < 
to < 1, PM covers the entire system and there exists no 
nematic state. Finally, if m = 1, we obtain FM again. 
These results are shown in FigfS] (b) . 

We remark here that the phase diagrams in the inset 
of Fig. |4]and Fig. [5] are analogous to the superfiuid case 
given in Ref.[4j. In mean field theories both the lattice 
and superfiuid cases yield mean field energy of the same 
forms due to symmetry. 



we have phase separation since the phase transition is 
first-order. If x is the fraction of nematic state, then 
x = 1 — m. This phase diagram is drawn in FigJS] (a). 

(b) < 4Q < Fig 0] shows AE for different states. 
AE(8i) is shown with a red dashed line, while AE(9 ) 
with a blue line as in Fig. [3J We can see that the red 
dashed line crosses zero at A n , p defined as (|53|) and then 
merges to the blue line at the point 



A 



pj 



zA. 



(56) 



In the regime: < A < A„ !P the ground state is nematic, 
for X n p < A < X p f the system is partially magnetic and 
one has a fully magnetic state if A > A P) /. Therefore 
the system undergoes two phase transitions: a first-order 
phase transition from NM to PM at X n>p and a second- 
order phase transition from PM to FM at A p ./. The 
second phase transition is second order due to the fact 
that 6*i goes to zero while A approaching zA, and there- 
fore the transition is continuous for the order parameter. 
This yields the phase diagram in the inset of Fig. [H 

The same question arises: if we have a net averaged 
magnetization per site to, which state we will achieve. 
To see this, we have to calculate the net magnetization 
for the PM state. From Eq.(T47|, (|48]) and gU) we can 
calculate to 



which leads to 



I |2 ^2 

171 =T 



TO = = COS Wi = — — 

h zA 



(57) 



(58) 



by using Eq. (|43|) . Therefore in the regime of PM states, 
A n ,p < A < Ap./, the magnetization lies in the region 




< to < 1. 



The PM state has a minimum magnetization 



(59) 



(60) 



As a result, if m = 
For < to < m m i n 



0, the system is purely nematic. 
phase separation occurs. One has 



B. Finite Temperature 

Before we determine the phase diagram for finite tem- 
perature, we first figure out different phases by investi- 
gating eigenstates of the density matrix p (fTTj) . After 
diagonalizing it, p is in its diagonalized form 

p = P + |*+)(* + | + P_|tf_)<*_| + P z \z)(z\, (61) 

where 



(62) 



p± = po ± y pI + pi. P z = Pf 

The eigenvectors read 

|*±)=u±|a:)+»±|tf), (63) 

where 

D 



W_ = V-L = 



ip2 

Vd^TpJ 



(64) 



with the definition: 



d = p3 + \ pj + pi 



(65) 



In order to obtain the true phases we have to solve 
the self-consistent equations Eqs. (f2"5|) to Eqs. (|3"0|) . As 
discussed in the last section, one can categorize these 
self-consistent equations into two groups: (1) /13 = and 
(2) /13 5^ with a constant h%. In the case /13 = 0, 
Pxx = Pyyii-e.hz — 0), the eigenvalues (|62p reads 



P± = P0 ± Pi = Pxx ± \Pxy\- 



(66) 



According to Eq. 



|^+) thus has the form 



u + = V- = ^7= and u_ = v + = 



|* + ) = — (Ix) + %)) = -| t), 



while | "F^) reads 



V2 



(\x)-i\y))=i\l). 



(67) 



(68) 
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— T/zlKI=0.1 
■■■■ T/zlKI=0.15 

- - T/zlKI = 0.2 
■-■ T/zlKI = 0.25 



m 2« 
0.8 - 

m m,n<- 
0.6 - 



m 



T/zlKI = 0.l 
T/zlKI = 0.15 
T/zlKI = 0.2 
T/zlKI = 0.25 



FIG. 6: (color online)The average magnetization per site m 
vs. A for different temperatures The parameters: Q/z\K\ — 
0.1, J/K = 0.91, A = §.\\K\. We can see that there is a huge 
magnetization jump at A„/(~ 0.446) even for T/z\K\ = 0.25. 
The maximum magnetization in NM is mi and the minimum 
magnetization in FM is m,2. 



FIG. 7: (color online) The average magnetization per site m 
vs. A for different temperatures. Q/z\K\ = 0.005 and other 
parameters are same as Fig. [6] m m i„ is the minimum mag- 
netization at which PM exists. A first-order phase transition 
occurs at X n , p ~ 0.06 from NM to PM and a second-order 
phase transition occurs at X p j depending on the tempera- 
ture. 



Therefore the system is a mixed state of | f), | I) and \z). 

In the second case that /13 7^ and /12 = const., the 
eigenstate \^+) can be rewritten as 

|*+) = a|T)+7U) (69) 

where 

a = --j=(u + -iv + ), 7= -^=(u+ (70) 

a, 7 are nonzero real numbers and 0^7. Similarly, | \I/ ) 

has the form 

l*-)=-7lt>+a|4>- (71) 

We can easily prove that |^+) is orthonormal to 

The system is a mixed state with (a| f ) + 7I J,)), (— 7I f 

) +a\ J.)) and \z). 

Numerically we solved the self-consistent equations 
and calculated their free energy according to Eqs. (|3"Tj) . 
dSU), (O and d35). J/K = 0.91 has been used to be 
close to those for 23 Na. In this case A = 0.091|it"j. In 
order to find the convergent solution quickly, we start 
with low temperature (T / z\K\ = 0.1) and extends the 
temperature step by step by using the final results as an 
initial input for the next temperature. We have evalu- 
ated the phase diagram up to T / z\K\ = 0.25. We note 
that if A = Q = 0, the nematic state becomes disordered 
at T / z\K\ ~ 0.36. We illustrate our result with two Q 
values: Q/z\K\ = 0.1 and Q/z\K\ = 0.005 to represent 
two regimes as for Fig. [3] and Fig. [4j We separate the 
two sets of self-consistent solutions: one with /13 = and 
one with /13 =/= and hi — const.. For the set of zero /13, 
two subsets occur. The first one contains the points with 
small fi2'. h% -C 1 and hi of the second subset is of or- 
der 1. Compared with the solutions of zero temperature, 
the first subset is a continuous evolution with the tem- 
perature T from the nematic solution, therefore we can 



still call these solution nematic (NM) , while in of case of 
large hi the solutions correspond to the fully magnetic 
states (FM) at zero temperature. On the other hand, if 
/i3 7^ and I12 — const., the states we obtain evolve from 
the PM states at zero temperature, we can still call them 
partially magnetic. 

(a) Q/z\K\ = 0.1 : in this case, we can calculate the 
free energy vs. A for the three different sets discussed 
above. The result is very similar to Fig. [3] for each tem- 
perature except that the free energy for nematic phase is 
not constant anymore but a monotonic decreasing func- 
tion of A. The transition points X n j of the first order 
phase transition stay almost the same for all tempera- 
tures, for this Q, X n j ~ 0.446. The nematic phase at fi- 
nite temperatures is not a pure state anymore, it contains 
mostly the nematic state \z) and with small amounts of 
I X) and I I) due to the fact that 7^ /12 "C 1. Therefore 
the magnetization is not zero. On the other hand, the 
FM phase is a mixture of large amount of | j) and small 
amounts of | [) and \z), as a result, the magnetization is 
smaller than one. PM phase can not appear here. 

We can also calculate magnetization for NM and FM 
phases by the relation 

m = Tr(n T - n^p. (72) 

It yields m = 2p 2 - Fig- El shows the magnetization at 
different temperatures in terms of A. to increases mono- 
tonically till it reaches its maximum mi at X n j and then 
jumps to value mi- At the end it increases to the fully 
magnetic state m = 1 if A ^> 1. At higher temperatures, 
77ii increases and mi decreases due to the fact that NM 
and FM mix more and more different states. As a result, 
if < to < toi, the system can be a uniform NM phase, 
while TO2 < J7i < 1, we obtain a uniform FM phase. In 
between, mi < to < to 2 , NM and FM coexist and they 
are phase separated since the phase transition is first- 



(a) 



NM 



NMandFM , P S 



FM 



m. 



m, 1 



NM NMandPM , P.S. PM FM 

(b) * • 



a m, 



FIG. 8: Phase diagrams for finite temperatures in terms of 
average net magnetization m. (a) large Q (b) small Q. P.S. 
means phase separation. This figure corresponds to Fig. [3] 
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order. We summarize the result in Fig. [5J (a). 

(b) Q/z\K\ = 0.005: similar to the zero tempera- 
ture case, PM phase appears here. The free energy 
curves vs. A for different temperatures are similar to 
Fig. [4j The first-order phase transition point stays 
the same: X np ~ 0.06 for all temperatures, while the 
second-order phase-transition point X p f changes: X p f — 
0.091,0.09,0.088,0.084 for T/z\K\ = 0.1,0.15,0.2,0.25, 
respectively. The reason is that for PM phase the con- 
straint hi < h has to be satisfied, it demands 



A < zA- 



2smh(3h 



2 cosh (ih 



K ,f(T). 



(73) 



At zero temperature, A Pj /(0) = zA which agrees with 
the result we obtained in the last section. With increas- 
ing T (decreasing /3), e^ h '° is getting larger, A p / is thus 
decreasing. 

As discussed above, PM phase has large amount of 
(a | t) +7l I)) mixed with small amounts of (— 7I f ) +a\ j 
) ) and I z) . The magnetization of the system for differ- 
ent temperatures as a function of A is shown in Fig. [7] 
In NM phase, m increases and then jumps to m m i n at 



A,, 



The system undergoes a first order phase transi- 



tion. In PM phase m ascends to 777,2 and the system 
changes continuously to FM phase. We conclude that if 
1 < 777 < 777i a homogeneous NM phase is achievable in 
experiments. In the case that mi < m < m m i n , NM 
and PM phases coexist but separate spatially. In the 
regime: m m j„ < m < 7772, PM phase with different mag- 
netization is the only phase in the system. In the end, if 
7772 < 777 < 1, we obtain FM phase. The phase diagram 
is plotted in Fig. [51(b). 



V. DISCUSSION AND CONCLUSION 

In the case A = and Q = 0, there exists a first- 
order phase transition between the nematic state and the 
disordered state at T/z\K\ = 0.36 [22]. The question 
arises naturally that if Q is nonzero, how the first-order 
phase transition develops. Fig. [5] shows the density p z 
as a function of temperature with increasing Q. Note 
that p z = 1/3 corresponds to a state with 0(3) sym- 
metry. We found that the first-order phase transition 
exists till Q = 0.0022:1^1 (red-dotted line) and then it 



FIG. 9: (color online)p z vs. temperature for different Q's. In 
this figure, A = 0. The first-order phase transition predicted 
for Q = still exists till Q = Q.QQ2z\K\. Even for larger Q, 
we can still see a rapid decrease of p z when T increases from 
the low to high temperature regime. 



turns to be a sharp crossover even till Q ~ 0.02z|A"| 
(thick yellow dashed line). From Ref.Q, Q is related 
to the magnetic field B: 4Q = 278 * B 2 (HzG~ 2 ), that 
means for Q = 0.002z|iY|, B = 0.0134^1^0 with 
Hz as the unit of \K\ . If the superexchange parame- 
ter is \K\ ~ lOOffzfil, this yields B = 0.134G. For 
Q = 0.02z\K\, B = 0.42G. Experimentally one can 
reach B < 0.01G, therefore the first-order phase tran- 
sition and the sharp crossover can be observable. On the 
other hand, as shown in text, the first-order phase transi- 
tion from the NM state to FM state for large Q and from 
the NM state to PM state for smaller Q remain at fi- 
nite temperature. We conclude that these phases would 
phase separate into different spatial regions. One may 
ask whether, instead of phase separation, one can have, 
for example, the ferromagnetic sites appear in the form 
of linear or planar stripes within the nematic regions. We 
exclude this for the following reason. According to Eq. 
@ and U2 — Uq for 23 Na, we obtain e ~ e 2 < e\ = 0, 
where e\ is the energy for total spin for two atoms equal 
to 1. Consider two neighboring sites. From the Clebsch- 
Gordan coefficients we can write down |0,0) as a linear 
combination of \F to t = 0) and \F to t — 2), and |1, 1) only 
exists in \F to t — 2), while |1,0) and |0, 1) must involve 
the high energy \F to t = 1) state. Therefore it costs more 
energy if the system builds a domain wall than just put 
the same m/ state as neighbors. That is the reason why 
the system prefers a spatially separated phase than stripe 
phases. A stripe phase is not favored because it needs to 
build more than one domain wall. 

To conclude, we have shown the phase diagrams for 
a spin-1 polar Bose gas loaded in an strongly repulsive 
optical lattice. There exist three different phases: the ne- 
matic (NM), fully magnetic (PM) and partially magnetic 
(PM) phases depending on the parameter regime of the 
system. A first-order phase transition from NM to FM 
or from NM to PM has been predicted. A second-order 
phase transition from PM to FM is also found. These 
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phase transitions are robust even at finite temperatures. 
Therefore they should be observable in experiments. 



APPENDIX A: DENSITY MATRIX 



To obtain p defined as Eq. l|22|) . we have to calculate 
e -P H eff first, ft can be written in the form 



,-0H = e /3(E?=o hm+hzT,) 



(Al) 



Since t z commutes with all <Xi, we can take e l3h ^ a ^ out of 
the exponential. One can prove the relation with prop- 
erties of Pauli matrices: 



(A2) 



where h is a three dimensional normal vector and a ~ 
[ci, (72, 03] T . By using this, Eq. (|Alj) has the form 



e -0H =e 0h, Tz + e 0ho CQsh p h aQ 



T a 2 + -a 3 , e 



0h o 



sinh (3h. 



(A3) 



It yields 

Tr e' 13 " = e ph * + 2e ph ° cosh (3h. (A4) 
Eq. ([13]) together with Eq. ([M]) are thus obtained. 

APPENDIX B: FREE ENERGY 



The internal energy can be calculated by using Eq. {T 
as follows 

1 



E int =-TrpH° eff 

= \z (K P aa + f)p° 



(Bl) 



+ z [J 



K 



{p " p u + p p • j 



This term can be simplified to a form 



zK 
~2~ 



E mt = ^ Trp 2 + 2z{J~K)p 2 2 + 



(B2) 



Note that Tr p 2 = 2p 2 + 2^2 + 2^2 + p 2 . Inserting Eq. J24]) 
into Eq. (|B2|) . we obtain internal energy as Eq. (|32jl . 
-Beat can be calculated in a similar way: 



E ext =Trp{4Q(\x){x\ + \y){y\) 
+iX(\x)(y\-\y)(x\)} 
=4Q - AQp zz - 2i\p xv , 



(B3) 



this yields 



E ext = 4Q - 4Qp z - 2\p 2 . (B4) 
Therefore Eq. ([33|) is obtained by using Eq. (124]) . 



Finally — TS* can be obtained by the definition of en- 
tropy: 



S — —kB Trpln p. 



(B5) 



It yields 



-TS = - Tr pH, 



eff 



ilnTre^ 

P 



--(2p - l)h + 2p 2 h 2 + 2p 3 h 3 + p z h z (B6) 
_ ln (e/">« +2cosh/3/i), 



and therefore Eq. 
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